

#include "slepbase.h"

void experiment_submatmult(){
	slep_matrix *A,*ATy;
	slep_vector *y;
	int i;
	int ind[] = {0,4,5,7};
	int k=3;

	y = slepAllocDenseVector(7);
	A = slepAllocDenseMatrix(7,9);
	ATy = slepAllocDenseMatrix(9,3);

	DMEMSET(y->data.vctr,1.0,7,i);
	DMEMSET(A->data.mtrx,1.0,7*9,i);

//	slep_randn(y->data.vctr,7,0.0,1.0);
//	slep_randn(A->data.mtrx,7*9,0.0,1.0);

	DMEMSET(ATy->data.mtrx,0.0,9*3,i);
	slepVerbConsoleWrDMat(ATy);	
	{
		slep_matrix A_part = *A;
		slep_vector ATy_i ;
		slep_vector y_i;
		ATy_i.dim = A->n; 
		ATy_i.inc = 1;
		ATy_i.sparsity = SlepDense;
		y_i.inc = 1;
		y_i.sparsity = SlepDense;
		for(i=0; i<k; i++){
			A_part.data.mtrx = &(A->data.mtrx[ind[i]]);
			A_part.m = ind[i+1] - ind[i];
			ATy_i.data.vctr = &(ATy->data.mtrx[i*ATy->ld]);
			y_i.data.vctr = &(y->data.vctr[ind[i]]);
			y_i.dim = ind[i+1] - ind[i];
			slep_daAxpby(SlepTranspose,1.0,&A_part, &y_i, 0.0, &ATy_i);
		}
	}
	slepVerbConsoleWrDMat(A);
	slepVerbConsoleWrDVec(y);
	slepVerbConsoleWrDMat(ATy);

}